home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_K0.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.6 KB  |  215 lines
  1. /* specfunc/bessel_K0.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_bessel.h>
  27.  
  28. #include "error.h"
  29.  
  30. #include "chebyshev.h"
  31. #include "cheb_eval.c"
  32.  
  33. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34.  
  35. /* based on SLATEC bk0(), bk0e() */
  36.  
  37. /* chebyshev expansions 
  38.  
  39.  series for bk0        on the interval  0.        to  4.00000d+00
  40.                     with weighted error   3.57e-19
  41.                      log weighted error  18.45
  42.                    significant figures required  17.99
  43.                     decimal places required  18.97
  44.  
  45.  series for ak0        on the interval  1.25000d-01 to  5.00000d-01
  46.                     with weighted error   5.34e-17
  47.                      log weighted error  16.27
  48.                    significant figures required  14.92
  49.                     decimal places required  16.89
  50.  
  51.  series for ak02       on the interval  0.        to  1.25000d-01
  52.                     with weighted error   2.34e-17
  53.                      log weighted error  16.63
  54.                    significant figures required  14.67
  55.                     decimal places required  17.20
  56. */
  57.  
  58. static double bk0_data[11] = {
  59.   -0.03532739323390276872,
  60.    0.3442898999246284869, 
  61.    0.03597993651536150163,
  62.    0.00126461541144692592,
  63.    0.00002286212103119451,
  64.    0.00000025347910790261,
  65.    0.00000000190451637722,
  66.    0.00000000001034969525,
  67.    0.00000000000004259816,
  68.    0.00000000000000013744,
  69.    0.00000000000000000035
  70. };
  71. static cheb_series bk0_cs = {
  72.   bk0_data,
  73.   10,
  74.   -1, 1,
  75.   10
  76. };
  77.  
  78. static double ak0_data[17] = {
  79.   -0.07643947903327941,
  80.   -0.02235652605699819,
  81.    0.00077341811546938,
  82.   -0.00004281006688886,
  83.    0.00000308170017386,
  84.   -0.00000026393672220,
  85.    0.00000002563713036,
  86.   -0.00000000274270554,
  87.    0.00000000031694296,
  88.   -0.00000000003902353,
  89.    0.00000000000506804,
  90.   -0.00000000000068895,
  91.    0.00000000000009744,
  92.   -0.00000000000001427,
  93.    0.00000000000000215,
  94.   -0.00000000000000033,
  95.    0.00000000000000005
  96. };
  97. static cheb_series ak0_cs = {
  98.   ak0_data,
  99.   16,
  100.   -1, 1,
  101.   10
  102. };
  103.  
  104. static double ak02_data[14] = {
  105.   -0.01201869826307592,
  106.   -0.00917485269102569,
  107.    0.00014445509317750,
  108.   -0.00000401361417543,
  109.    0.00000015678318108,
  110.   -0.00000000777011043,
  111.    0.00000000046111825,
  112.   -0.00000000003158592,
  113.    0.00000000000243501,
  114.   -0.00000000000020743,
  115.    0.00000000000001925,
  116.   -0.00000000000000192,
  117.    0.00000000000000020,
  118.   -0.00000000000000002
  119. };
  120. static cheb_series ak02_cs = {
  121.   ak02_data,
  122.   13,
  123.   -1, 1,
  124.   8
  125. };
  126.  
  127.  
  128. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  129.  
  130. int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result)
  131. {
  132.   /* CHECK_POINTER(result) */
  133.  
  134.   if(x <= 0.0) {
  135.     DOMAIN_ERROR(result);
  136.   }
  137.   else if(x <= 2.0) {
  138.     const double lx = log(x);
  139.     const double ex = exp(x);
  140.     int stat_I0;
  141.     gsl_sf_result I0;
  142.     gsl_sf_result c;
  143.     cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);
  144.     stat_I0 = gsl_sf_bessel_I0_e(x, &I0);
  145.     result->val  = ex * ((-lx+M_LN2)*I0.val - 0.25 + c.val);
  146.     result->err  = ex * ((M_LN2+fabs(lx))*I0.err + c.err);
  147.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  148.     return stat_I0;
  149.   }
  150.   else if(x <= 8.0) {
  151.     const double sx = sqrt(x);
  152.     gsl_sf_result c;
  153.     cheb_eval_e(&ak0_cs, (16.0/x-5.0)/3.0, &c);
  154.     result->val  = (1.25 + c.val) / sx;
  155.     result->err  = c.err / sx;
  156.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  157.     return GSL_SUCCESS;
  158.   }
  159.   else {
  160.     const double sx = sqrt(x);
  161.     gsl_sf_result c;
  162.     cheb_eval_e(&ak02_cs, 16.0/x-1.0, &c);
  163.     result->val  = (1.25 + c.val) / sx;
  164.     result->err  = (c.err + GSL_DBL_EPSILON) / sx;
  165.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  166.     return GSL_SUCCESS;
  167.   } 
  168. }
  169.  
  170.  
  171. int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result)
  172. {
  173.   /* CHECK_POINTER(result) */
  174.  
  175.   if(x <= 0.0) {
  176.     DOMAIN_ERROR(result);
  177.   }
  178.   else if(x <= 2.0) {
  179.     const double lx = log(x);
  180.     int stat_I0;
  181.     gsl_sf_result I0;
  182.     gsl_sf_result c;
  183.     cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);
  184.     stat_I0 = gsl_sf_bessel_I0_e(x, &I0);
  185.     result->val  = (-lx+M_LN2)*I0.val - 0.25 + c.val;
  186.     result->err  = (fabs(lx) + M_LN2) * I0.err + c.err;
  187.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  188.     return stat_I0;
  189.   }
  190.   else {
  191.     gsl_sf_result K0_scaled;
  192.     int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled);
  193.     int stat_e  = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x),
  194.                                            K0_scaled.val, K0_scaled.err,
  195.                        result);
  196.     return GSL_ERROR_SELECT_2(stat_e, stat_K0);
  197.   }
  198. }
  199.  
  200.  
  201. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  202.  
  203. #include "eval.h"
  204.  
  205. double gsl_sf_bessel_K0_scaled(const double x)
  206. {
  207.   EVAL_RESULT(gsl_sf_bessel_K0_scaled_e(x, &result));
  208. }
  209.  
  210. double gsl_sf_bessel_K0(const double x)
  211. {
  212.   EVAL_RESULT(gsl_sf_bessel_K0_e(x, &result));
  213. }
  214.  
  215.